# Loading the Drive helper and mount your Google Drive as a drive in the virtual machine
from google.colab import drive
drive.mount('/content/drive')
Mounted at /content/drive
# Installing some libraries that are not on Colab by default
!pip install rasterio
!pip install geopandas
!pip install rasterstats
!pip install earthengine-api
!pip install requests
!pip install sentinelsat
# Importing libraries
import geopandas as gpd
import rasterio
from rasterio import plot
from rasterio.plot import show_hist
import matplotlib.pyplot as plt
import numpy as np
from osgeo import gdal, ogr, osr
import json
import os
from os import listdir
from os.path import isfile, isdir, join
import math
from pprint import pprint
import shutil
import sys
import zipfile
import requests
import io
import webbrowser
import ee
import pandas as pd
import rasterio.mask
import xarray as xr
import matplotlib.colors as colors
import matplotlib
import glob
from rasterio.plot import show
from rasterio.mask import mask
from shapely.geometry import mapping
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# defining the root directory where our data are to be stored
rootdir = '/content/drive/MyDrive/satelliteCW1' # this is where pygge.py is saved on my Google Drive
if rootdir not in sys.path:
sys.path.append(rootdir)
# importing the pygge library of functions for this module
import pygge
%matplotlib inline
Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Collecting rasterio
Downloading rasterio-1.3.6-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (20.0 MB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 20.0/20.0 MB 104.8 MB/s eta 0:00:00
Collecting affine (from rasterio)
Downloading affine-2.4.0-py3-none-any.whl (15 kB)
Requirement already satisfied: attrs in /usr/local/lib/python3.10/dist-packages (from rasterio) (23.1.0)
Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from rasterio) (2022.12.7)
Requirement already satisfied: click>=4.0 in /usr/local/lib/python3.10/dist-packages (from rasterio) (8.1.3)
Collecting cligj>=0.5 (from rasterio)
Downloading cligj-0.7.2-py3-none-any.whl (7.1 kB)
Requirement already satisfied: numpy>=1.18 in /usr/local/lib/python3.10/dist-packages (from rasterio) (1.22.4)
Collecting snuggs>=1.4.1 (from rasterio)
Downloading snuggs-1.4.7-py3-none-any.whl (5.4 kB)
Collecting click-plugins (from rasterio)
Downloading click_plugins-1.1.1-py2.py3-none-any.whl (7.5 kB)
Requirement already satisfied: setuptools in /usr/local/lib/python3.10/dist-packages (from rasterio) (67.7.2)
Requirement already satisfied: pyparsing>=2.1.6 in /usr/local/lib/python3.10/dist-packages (from snuggs>=1.4.1->rasterio) (3.0.9)
Installing collected packages: snuggs, cligj, click-plugins, affine, rasterio
Successfully installed affine-2.4.0 click-plugins-1.1.1 cligj-0.7.2 rasterio-1.3.6 snuggs-1.4.7
Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Collecting geopandas
Downloading geopandas-0.13.0-py3-none-any.whl (1.1 MB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 1.1/1.1 MB 18.6 MB/s eta 0:00:00
Collecting fiona>=1.8.19 (from geopandas)
Downloading Fiona-1.9.4-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (16.5 MB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 16.5/16.5 MB 36.6 MB/s eta 0:00:00
Requirement already satisfied: packaging in /usr/local/lib/python3.10/dist-packages (from geopandas) (23.1)
Requirement already satisfied: pandas>=1.1.0 in /usr/local/lib/python3.10/dist-packages (from geopandas) (1.5.3)
Collecting pyproj>=3.0.1 (from geopandas)
Downloading pyproj-3.5.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (7.7 MB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 7.7/7.7 MB 124.7 MB/s eta 0:00:00
Requirement already satisfied: shapely>=1.7.1 in /usr/local/lib/python3.10/dist-packages (from geopandas) (2.0.1)
Requirement already satisfied: attrs>=19.2.0 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (23.1.0)
Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (2022.12.7)
Requirement already satisfied: click~=8.0 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (8.1.3)
Requirement already satisfied: click-plugins>=1.0 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (1.1.1)
Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (0.7.2)
Requirement already satisfied: six in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (1.16.0)
Requirement already satisfied: python-dateutil>=2.8.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.0->geopandas) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.0->geopandas) (2022.7.1)
Requirement already satisfied: numpy>=1.21.0 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.0->geopandas) (1.22.4)
Installing collected packages: pyproj, fiona, geopandas
Successfully installed fiona-1.9.4 geopandas-0.13.0 pyproj-3.5.0
Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Collecting rasterstats
Downloading rasterstats-0.18.0-py3-none-any.whl (17 kB)
Requirement already satisfied: affine<3.0 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (2.4.0)
Requirement already satisfied: click>7.1 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (8.1.3)
Requirement already satisfied: cligj>=0.4 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (0.7.2)
Collecting fiona<1.9 (from rasterstats)
Downloading Fiona-1.8.22-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (16.6 MB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 16.6/16.6 MB 117.5 MB/s eta 0:00:00
Requirement already satisfied: numpy>=1.9 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (1.22.4)
Requirement already satisfied: rasterio>=1.0 in /usr/local/lib/python3.10/dist-packages (from rasterstats) (1.3.6)
Collecting simplejson (from rasterstats)
Downloading simplejson-3.19.1-cp310-cp310-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl (137 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 137.9/137.9 kB 22.3 MB/s eta 0:00:00
Requirement already satisfied: shapely in /usr/local/lib/python3.10/dist-packages (from rasterstats) (2.0.1)
Requirement already satisfied: attrs>=17 in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (23.1.0)
Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (2022.12.7)
Requirement already satisfied: click-plugins>=1.0 in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (1.1.1)
Requirement already satisfied: six>=1.7 in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (1.16.0)
Collecting munch (from fiona<1.9->rasterstats)
Downloading munch-3.0.0-py2.py3-none-any.whl (10 kB)
Requirement already satisfied: setuptools in /usr/local/lib/python3.10/dist-packages (from fiona<1.9->rasterstats) (67.7.2)
Requirement already satisfied: snuggs>=1.4.1 in /usr/local/lib/python3.10/dist-packages (from rasterio>=1.0->rasterstats) (1.4.7)
Requirement already satisfied: pyparsing>=2.1.6 in /usr/local/lib/python3.10/dist-packages (from snuggs>=1.4.1->rasterio>=1.0->rasterstats) (3.0.9)
Installing collected packages: simplejson, munch, fiona, rasterstats
Attempting uninstall: fiona
Found existing installation: Fiona 1.9.4
Uninstalling Fiona-1.9.4:
Successfully uninstalled Fiona-1.9.4
Successfully installed fiona-1.8.22 munch-3.0.0 rasterstats-0.18.0 simplejson-3.19.1
Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Requirement already satisfied: earthengine-api in /usr/local/lib/python3.10/dist-packages (0.1.350)
Requirement already satisfied: google-cloud-storage in /usr/local/lib/python3.10/dist-packages (from earthengine-api) (2.8.0)
Requirement already satisfied: google-api-python-client>=1.12.1 in /usr/local/lib/python3.10/dist-packages (from earthengine-api) (2.84.0)
Requirement already satisfied: google-auth>=1.4.1 in /usr/local/lib/python3.10/dist-packages (from earthengine-api) (2.17.3)
Requirement already satisfied: google-auth-httplib2>=0.0.3 in /usr/local/lib/python3.10/dist-packages (from earthengine-api) (0.1.0)
Requirement already satisfied: httplib2<1dev,>=0.9.2 in /usr/local/lib/python3.10/dist-packages (from earthengine-api) (0.21.0)
Requirement already satisfied: requests in /usr/local/lib/python3.10/dist-packages (from earthengine-api) (2.27.1)
Requirement already satisfied: google-api-core!=2.0.*,!=2.1.*,!=2.2.*,!=2.3.0,<3.0.0dev,>=1.31.5 in /usr/local/lib/python3.10/dist-packages (from google-api-python-client>=1.12.1->earthengine-api) (2.11.0)
Requirement already satisfied: uritemplate<5,>=3.0.1 in /usr/local/lib/python3.10/dist-packages (from google-api-python-client>=1.12.1->earthengine-api) (4.1.1)
Requirement already satisfied: cachetools<6.0,>=2.0.0 in /usr/local/lib/python3.10/dist-packages (from google-auth>=1.4.1->earthengine-api) (5.3.0)
Requirement already satisfied: pyasn1-modules>=0.2.1 in /usr/local/lib/python3.10/dist-packages (from google-auth>=1.4.1->earthengine-api) (0.3.0)
Requirement already satisfied: six>=1.9.0 in /usr/local/lib/python3.10/dist-packages (from google-auth>=1.4.1->earthengine-api) (1.16.0)
Requirement already satisfied: rsa<5,>=3.1.4 in /usr/local/lib/python3.10/dist-packages (from google-auth>=1.4.1->earthengine-api) (4.9)
Requirement already satisfied: pyparsing!=3.0.0,!=3.0.1,!=3.0.2,!=3.0.3,<4,>=2.4.2 in /usr/local/lib/python3.10/dist-packages (from httplib2<1dev,>=0.9.2->earthengine-api) (3.0.9)
Requirement already satisfied: google-cloud-core<3.0dev,>=2.3.0 in /usr/local/lib/python3.10/dist-packages (from google-cloud-storage->earthengine-api) (2.3.2)
Requirement already satisfied: google-resumable-media>=2.3.2 in /usr/local/lib/python3.10/dist-packages (from google-cloud-storage->earthengine-api) (2.5.0)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in /usr/local/lib/python3.10/dist-packages (from requests->earthengine-api) (1.26.15)
Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.10/dist-packages (from requests->earthengine-api) (2022.12.7)
Requirement already satisfied: charset-normalizer~=2.0.0 in /usr/local/lib/python3.10/dist-packages (from requests->earthengine-api) (2.0.12)
Requirement already satisfied: idna<4,>=2.5 in /usr/local/lib/python3.10/dist-packages (from requests->earthengine-api) (3.4)
Requirement already satisfied: googleapis-common-protos<2.0dev,>=1.56.2 in /usr/local/lib/python3.10/dist-packages (from google-api-core!=2.0.*,!=2.1.*,!=2.2.*,!=2.3.0,<3.0.0dev,>=1.31.5->google-api-python-client>=1.12.1->earthengine-api) (1.59.0)
Requirement already satisfied: protobuf!=3.20.0,!=3.20.1,!=4.21.0,!=4.21.1,!=4.21.2,!=4.21.3,!=4.21.4,!=4.21.5,<5.0.0dev,>=3.19.5 in /usr/local/lib/python3.10/dist-packages (from google-api-core!=2.0.*,!=2.1.*,!=2.2.*,!=2.3.0,<3.0.0dev,>=1.31.5->google-api-python-client>=1.12.1->earthengine-api) (3.20.3)
Requirement already satisfied: google-crc32c<2.0dev,>=1.0 in /usr/local/lib/python3.10/dist-packages (from google-resumable-media>=2.3.2->google-cloud-storage->earthengine-api) (1.5.0)
Requirement already satisfied: pyasn1<0.6.0,>=0.4.6 in /usr/local/lib/python3.10/dist-packages (from pyasn1-modules>=0.2.1->google-auth>=1.4.1->earthengine-api) (0.5.0)
Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Requirement already satisfied: requests in /usr/local/lib/python3.10/dist-packages (2.27.1)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in /usr/local/lib/python3.10/dist-packages (from requests) (1.26.15)
Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.10/dist-packages (from requests) (2022.12.7)
Requirement already satisfied: charset-normalizer~=2.0.0 in /usr/local/lib/python3.10/dist-packages (from requests) (2.0.12)
Requirement already satisfied: idna<4,>=2.5 in /usr/local/lib/python3.10/dist-packages (from requests) (3.4)
Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Collecting sentinelsat
Downloading sentinelsat-1.2.1-py3-none-any.whl (48 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 48.8/48.8 kB 3.2 MB/s eta 0:00:00
Requirement already satisfied: requests in /usr/local/lib/python3.10/dist-packages (from sentinelsat) (2.27.1)
Requirement already satisfied: click>=7.1 in /usr/local/lib/python3.10/dist-packages (from sentinelsat) (8.1.3)
Collecting html2text (from sentinelsat)
Downloading html2text-2020.1.16-py3-none-any.whl (32 kB)
Collecting geojson>=2 (from sentinelsat)
Downloading geojson-3.0.1-py3-none-any.whl (15 kB)
Requirement already satisfied: tqdm>=4.58 in /usr/local/lib/python3.10/dist-packages (from sentinelsat) (4.65.0)
Collecting geomet (from sentinelsat)
Downloading geomet-1.0.0-py3-none-any.whl (28 kB)
Requirement already satisfied: six in /usr/local/lib/python3.10/dist-packages (from geomet->sentinelsat) (1.16.0)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in /usr/local/lib/python3.10/dist-packages (from requests->sentinelsat) (1.26.15)
Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.10/dist-packages (from requests->sentinelsat) (2022.12.7)
Requirement already satisfied: charset-normalizer~=2.0.0 in /usr/local/lib/python3.10/dist-packages (from requests->sentinelsat) (2.0.12)
Requirement already satisfied: idna<4,>=2.5 in /usr/local/lib/python3.10/dist-packages (from requests->sentinelsat) (3.4)
Installing collected packages: html2text, geomet, geojson, sentinelsat
Successfully installed geojson-3.0.1 geomet-1.0.0 html2text-2020.1.16 sentinelsat-1.2.1
A shapefile of Larache province is in my Google Drive. I drew a polygon and saved it as a shapefile on http://www.geojson.io.
# Connecting to Google Earth Engine API
# This opens a web page where i enter my account information and a verification code is also provided. This code is what i paste into the terminal.
!earthengine authenticate
ee.Initialize()
To authorize access needed by Earth Engine, open the following URL in a web browser and follow the instructions. If the web browser does not start automatically, please manually browse the URL below.
https://code.earthengine.google.com/client-auth?scopes=https%3A//www.googleapis.com/auth/earthengine%20https%3A//www.googleapis.com/auth/devstorage.full_control&request_id=m6jhmZYUMFpwsVc7HGY_G9JPMbQ8rSlu38E6KuTfon8&tc=qv5LaParODYGIYSxAjbHdPXmk0FuKvtCKGY7nEk4YsE&cc=It22cv2Oj5kPlxLjRdYc79Gu7DVH3ZkLsxjBNxXapOc
The authorization workflow will generate a code, which you should paste in the box below.
Enter verification code: 4/1AbUR2VM1jcMN2HebXQz4IYDnKP9_8Xg8RInpG8gTYo-ZrsMEsR5Pdi9gjSw
Successfully saved authorization token.
'''
--------------------------------------------------------------
The following block of code originates from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# set up your directories for the satellite data
# Note that we do all the downloading and data analysis on the temporary drive
# on Colab. We will copy the output directory to our Google Drive at the end.
# Colab has more disk space (about 40 GB free space) than Google Drive (15 GB).
# However, the data on the Colab disk space are NOT kept when you log out.
# path to your Google Drive
print("Connected to data directory: " + rootdir)
# path to your temporary drive on the Colab Virtual Machine. This will be removed
# when the session ends.
cd = "/content/work"
# directory for downloading the Sentinel-2 composites
# Note that we are using the 'join' function imported from the os library here
# It is an easy way of merging strings into a directory structure.
# It is clever and chooses the / or \ depending on whether you are on Windows or Linux.
downloaddir = join(cd, 'download') # where we save the downloaded images
# CAREFUL: This code removes the named directories and everything inside them to free up space
# Because the downloaddir is in your temporary drive, it should be empty at the
# start of the session.
# Note: shutil provides a lot of useful functions for file and directory management
try:
shutil.rmtree(downloaddir)
except:
print(downloaddir + " not found.")
# create the new directories, unless they already exist
os.makedirs(cd, exist_ok=True)
os.makedirs(downloaddir, exist_ok=True)
print("Connected to Colab temporary data directory: " + cd)
print("\nList of contents of " + rootdir)
for f in sorted(os.listdir(rootdir)):
print(f)
Connected to data directory: /content/drive/MyDrive/satelliteCW1 /content/work/download not found. Connected to Colab temporary data directory: /content/work List of contents of /content/drive/MyDrive/satelliteCW1 .ipynb_checkpoints 229010645_GY7709_CW1 BACKUP2.ipynb 229010645_GY7709_CW1 BACKUP3.ipynb 229010645_GY7709_CW1.html 229010645_GY7709_CW1.ipynb 229010645_GY7709_CW1_backup.ipynb 229010645_GY7709_CW1_original.html 229010645_GY7709_CW2.html 229010645_GY7709_CW2.ipynb NBR_larache_after.tiff NBR_larache_before.tiff NDVI_larache_after.tiff NDVI_larache_before.tiff __pycache__ backup burned_layers dNBR.tif dNBR.tiff dNBR_warped.tif larache_after larache_before larache_before_fire.tif larache_before_fire_warped.tif larache_new_shapefile larache_shapefiles leicestershire oakham practical_week32_SentinelSat.ipynb pygge.ipynb pygge.py sencredentials.txt taza unburned_layers
modifying some of the parameters and uploading the shapefile of Larache called "POLYGON.shp" .
'''
--------------------------------------------------------------
The following block of code is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# EDITING THE SEARCH OPTIONS BELOW
shapefile = join(rootdir, 'larache_new_shapefile', 'FULL_POLYGON.shp') # ESRI Shapefile of the study area
# checking whether the shapefile exists
if os.path.exists(shapefile):
print('Shapefile found: '+shapefile)
else:
print('ERROR: Shapefile not found: '+shapefile)
print('Upload a shapefile to your Google Drive directory: '+ rootdir)
# Defining a date range for the search
datefrom = '2022-07-01' # start date for imagery search
dateto = '2022-07-10' # end date for imagery search
time_range = [datefrom, dateto] # format as a list
# Defining which cloud cover percentage we accept in the images
clouds = 8 # maximum acceptable cloud cover in
Shapefile found: /content/drive/MyDrive/satelliteCW1/larache_new_shapefile/FULL_POLYGON.shp
'''
--------------------------------------------------------------
The following block of code originates from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# Getting the shapefile layer's extent, CRS and EPSG code
extent, outSpatialRef, epsg = pygge.get_shp_extent(shapefile)
print("Extent of the area of interest (shapefile):\n", extent)
print(type(extent))
print("\nCoordinate referencing system (CRS) of the shapefile:\n", outSpatialRef)
print('EPSG code: ', epsg)
Extent of the area of interest (shapefile):
(-6.08, -6.02, 35.21, 35.31)
<class 'tuple'>
Coordinate referencing system (CRS) of the shapefile:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AXIS["Latitude",NORTH],
AXIS["Longitude",EAST],
AUTHORITY["EPSG","4326"]]
EPSG code: 4326
Getting the extent of the shapefile into a format that Google Earth Engine understands.
'''
--------------------------------------------------------------
The following block of code originates from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# GEE needs a special format for defining an area of interest.
# It has to be a GeoJSON Polygon and the coordinates should be first defined in a list and then converted using ee.Geometry.
extent_list = list(extent)
print(extent_list)
print(type(extent_list))
# close the list of polygon coordinates by adding the starting node at the end again
# and make list elements in the form of coordinate pairs (y,x)
area_list = list([(extent[0], extent[2]),(extent[1], extent[2]),(extent[1], extent[3]),(extent[0], extent[3]),(extent[0], extent[2])])
print(area_list)
print(type(area_list))
search_area = ee.Geometry.Polygon(area_list)
print(search_area)
print(type(search_area))
[-6.08, -6.02, 35.21, 35.31]
<class 'list'>
[(-6.08, 35.21), (-6.02, 35.21), (-6.02, 35.31), (-6.08, 35.31), (-6.08, 35.21)]
<class 'list'>
ee.Geometry({
"functionInvocationValue": {
"functionName": "GeometryConstructors.Polygon",
"arguments": {
"coordinates": {
"constantValue": [
[
[
-6.08,
35.21
],
[
-6.02,
35.21
],
[
-6.02,
35.31
],
[
-6.08,
35.31
],
[
-6.08,
35.21
]
]
]
},
"evenOdd": {
"constantValue": true
}
}
}
})
<class 'ee.geometry.Geometry'>
Accessing the Sentinel-2 collection on Google Earth Engine and running our search.
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# Obtaining download links for image composites from an image collection on Google Earth Engine
# Name of the Sentinel 2 image collection
s2collection = ('COPERNICUS/S2')
# getting the median composite of Sentinel-2 images in the time range
s2median = pygge.obtain_image_sentinel(s2collection, time_range, search_area, clouds)
# Downloading the followig bands as a list of strings namely; Blue, Green,Red,NIR and SWIR
bands = [ 'B2','B3','B4', 'B8','B12']
print(bands)
# spatial resolution of the downloaded data
resolution = 20 # in units of metres
# Download images in Geotiff, using the get_url(name, image, scale, region) method
# ‘region’ is obtained from the area, but the format is adjusted using get_region(geom) method
search_region = pygge.get_region(search_area)
s2url = pygge.get_url('s2', s2median.select(bands), resolution, search_region, filePerBand=False)
print(s2url)
['B2', 'B3', 'B4', 'B8', 'B12'] https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/44e3d97d9adfca34655516a4c2fea849-11e933b4d0f41ebc8d8b877d0589bca7:getPixels
The next cell downloads the image composite as a zip file and unzips it.
'''
--------------------------------------------------------------
The following block of code originates from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# change directory to download directory
os.chdir(downloaddir)
# requesting information on the file to be downloaded
f = pygge.requests.get(s2url, stream =True)
# checking whether it is a zip file
check = zipfile.is_zipfile(io.BytesIO(f.content))
# either download the file as is, or unzip it
while not check:
f = requests.get(s2url, stream =True)
check = zipfile.is_zipfile(io.BytesIO(f.content))
else:
z = zipfile.ZipFile(io.BytesIO(f.content))
z.extractall()
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# where the downloaded Sentinel-2 images are stored
os.chdir(downloaddir)
print("contents of ", downloaddir, ":")
!ls -l
contents of /content/work/download : total 2120 -rw-r--r-- 1 root root 2168674 May 12 13:28 s2.tif
The downloaded file " s2.tif " is seen.
the downloaded images are saved to a temporary directory that will be deleted when the virtual machine is closed. To save the images to my local directory, this is how it went;
I went to my Google Colab folder in the panel on the left hand side.
Found the download directory and clicked on a Sentinel-2 image folder.
Right-clicked on it ,renamed it and selected 'download' to save it.
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# getting list of all tiff files in the directory
allfiles = [f for f in listdir(downloaddir) if isfile(join(downloaddir, f))]
print(allfiles)
# selecting the file for visualisation
thisfile = allfiles[0]
print(thisfile)
['s2.tif'] s2.tif
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# creating a figure
fig, (ax1) = plt.subplots(1,1, figsize=(7,7))
fig.patch.set_facecolor('white')
# the downloaded file is float32 data format
# for plotting, uint8 data format is needed
# plotting the image with full extent
pygge.easy_plot(thisfile, ax=ax1, bands=[3,2,1], percentiles=[0,99])
WARNING:rasterio._env:CPLE_AppDefined in s2.tif: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples. WARNING:rasterio._env:CPLE_AppDefined in TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# creating a figure with 2x3 subplots
fig, (ax1) = plt.subplots(1,1, figsize=(7,7))
fig.patch.set_facecolor('white')
# the downloaded file is float32 data format
# for plotting, uint8 data format is needed
# plotting the image with full extent
pygge.easy_plot(thisfile, ax=ax1, bands=[4,2,1], percentiles=[0,99])
WARNING:rasterio._env:CPLE_AppDefined in s2.tif: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples. WARNING:rasterio._env:CPLE_AppDefined in TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.
The coordinate reference system (CRS) of the downloaded image composite is not in the UK national map projection. hence it is reprojected.
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# printing the EPSG code of our shapefile into which we want to reproject the TCI images
print("Reprojecting image to EPSG projection ", epsg)
# making a file name for the new file
warpfile = thisfile.split(sep='.')[0] + '_warped.tif'
print("We are in this directory: ")
!pwd
print("Input file: ", thisfile)
print("Output file: ", warpfile)
# calling the easy_warp function
tmp = pygge.easy_warp(thisfile, warpfile, epsg)
Reprojecting image to EPSG projection 4326 We are in this directory: /content/work/download
WARNING:rasterio._env:CPLE_AppDefined in s2.tif: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples. WARNING:rasterio._env:CPLE_AppDefined in TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.
Input file: s2.tif Output file: s2_warped.tif Creating warped file:s2_warped.tif
To see the locations of our polygons on top of our image composite, we can do that with the Geopandas library.
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# creating a figure with subplots
fig, ax = plt.subplots(2,1, figsize=(10,16))
fig.patch.set_facecolor('white')
# plotting the image with full extent in true colour
pygge.easy_plot(warpfile, ax=ax[0], percentiles=[0,98], bands=[3,2,1],
shapefile=shapefile, fillcolor="none", linecolor="yellow",
title="LARACHE BEFORE FIRE")
# plotting the image with full extent in false colour
pygge.easy_plot(warpfile, ax=ax[1], percentiles=[0,98], bands=[4,2,1],
shapefile=shapefile, fillcolor="none", linecolor="yellow",
title="LARACHE BEFORE FIRE")
We change the search dates to get a new image after the fire had burnt out. Fires were first reported on 15 July 2022
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# EDITing THE SEARCH OPTIONS BELOW
shapefile = join(rootdir, 'larache_shapefiles', 'FULL_POLYGON.shp') # ESRI Shapefile of the study area
# checking whether the shapefile exists
if os.path.exists(shapefile):
print('Shapefile found: '+shapefile)
else:
print('ERROR: Shapefile not found: '+shapefile)
print('Upload a shapefile to your Google Drive directory: '+ rootdir)
# Defining a date range for the search
datefrom = '2022-07-30' # start date for imagery search
dateto = '2022-08-05' # end date for imagery search
time_range = [datefrom, dateto] # format as a list
# Defining percentage cloud cover accepted in the images
clouds = 10 # maximum acceptable cloud cover in %
Shapefile found: /content/drive/MyDrive/satelliteCW1/larache_shapefiles/FULL_POLYGON.shp
Get some information about our shapefile.
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# Getting the shapefile layer's extent, CRS and EPSG code
extent, outSpatialRef, epsg = pygge.get_shp_extent(shapefile)
print("Extent of the area of interest (shapefile):\n", extent)
print(type(extent))
print("\nCoordinate referencing system (CRS) of the shapefile:\n", outSpatialRef)
print('EPSG code: ', epsg)
Extent of the area of interest (shapefile):
(-6.08, -6.02, 35.21, 35.31)
<class 'tuple'>
Coordinate referencing system (CRS) of the shapefile:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AXIS["Latitude",NORTH],
AXIS["Longitude",EAST],
AUTHORITY["EPSG","4326"]]
EPSG code: 4326
Getting the extent of the shapefile into a format that Google Earth Engine understands.
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# GEE needs a special format for defining an area of interest.
# It has to be a GeoJSON Polygon and the coordinates should be first defined in a list and then converted using ee.Geometry.
extent_list = list(extent)
print(extent_list)
print(type(extent_list))
# close the list of polygon coordinates by adding the starting node at the end again
# and make list elements in the form of coordinate pairs (y,x)
area_list = list([(extent[0], extent[2]),(extent[1], extent[2]),(extent[1], extent[3]),(extent[0], extent[3]),(extent[0], extent[2])])
print(area_list)
print(type(area_list))
search_area = ee.Geometry.Polygon(area_list)
print(search_area)
print(type(search_area))
[-6.08, -6.02, 35.21, 35.31]
<class 'list'>
[(-6.08, 35.21), (-6.02, 35.21), (-6.02, 35.31), (-6.08, 35.31), (-6.08, 35.21)]
<class 'list'>
ee.Geometry({
"functionInvocationValue": {
"functionName": "GeometryConstructors.Polygon",
"arguments": {
"coordinates": {
"constantValue": [
[
[
-6.08,
35.21
],
[
-6.02,
35.21
],
[
-6.02,
35.31
],
[
-6.08,
35.31
],
[
-6.08,
35.21
]
]
]
},
"evenOdd": {
"constantValue": true
}
}
}
})
<class 'ee.geometry.Geometry'>
Accessing the Sentinel-2 collection on Google Earth Engine and running the search.
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# Obtain download links for image composites from an image collection on Google Earth Engine
# Name of the Sentinel 2 image collection
s2collection = ('COPERNICUS/S2')
# getting the median composite of Sentinel-2 images in the time range
s2median = pygge.obtain_image_sentinel(s2collection, time_range, search_area, clouds)
# Blue, Green,Red,NIR and SWIR are downloaded
bands = ['B2','B3','B4', 'B8','B12' ]
print(bands)
# spatial resolution of the downloaded data
resolution = 20 # in units of metres
# Downloading images in Geotiff, using the get_url(name, image, scale, region) method
# ‘region’ is obtained from the area, but the format has to be adjusted using get_region(geom) method
search_region = pygge.get_region(search_area)
s2url = pygge.get_url('s2', s2median.select(bands), resolution, search_region, filePerBand=False)
print(s2url)
['B2', 'B3', 'B4', 'B8', 'B12'] https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/3b62a644518dacef2cf989403d8c3e3f-44432a28d67b784f448895200db8318e:getPixels
The next cell downloads the image composite as a zip file and unzips it.
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# changing directory to download directory
os.chdir(downloaddir)
# requesting information on the file to be downloaded
f = pygge.requests.get(s2url, stream =True)
# checking whether it is a zip file
check = zipfile.is_zipfile(io.BytesIO(f.content))
# either download the file as is, or unzip it
while not check:
f = requests.get(s2url, stream =True)
check = zipfile.is_zipfile(io.BytesIO(f.content))
else:
z = zipfile.ZipFile(io.BytesIO(f.content))
z.extractall()
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# where we stored the downloaded Sentinel-2 images
os.chdir(downloaddir)
print("contents of ", downloaddir, ":")
!ls -l
contents of /content/work/download : total 2200 -rw-r--r-- 1 root root 2251399 Apr 27 11:10 s2.tif
the downloaded file is " s2.tif ".
the downloaded images are saved to a temporary directory that will be deleted when the virtual machine is closed. To save the images to my local directory, this is how it went;
I went to my Google Colab folder in the panel on the left hand side.
Found the download directory and clicked on a Sentinel-2 image folder.
Right-clicked on it ,renamed it and selected 'download' to save it.
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# getting list of all tiff files in the directory
allfiles = [f for f in listdir(downloaddir) if isfile(join(downloaddir, f))]
print(allfiles)
# selecting the file for visualisation
thisfile = allfiles[0]
print(thisfile)
['s2.tif'] s2.tif
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# creating a figure with 2x3 subplots
fig, (ax1) = plt.subplots(1,1, figsize=(7,7))
fig.patch.set_facecolor('white')
# the downloaded file is float32 data format
# for plotting, uint8 data format is needed
# plotting the image with full extent
pygge.easy_plot(thisfile, ax=ax1, bands=[3,2,1], percentiles=[0,99])
WARNING:rasterio._env:CPLE_AppDefined in s2.tif: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples. WARNING:rasterio._env:CPLE_AppDefined in TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# creating a figure with 2x3 subplots
fig, (ax1) = plt.subplots(1,1, figsize=(7,7))
fig.patch.set_facecolor('white')
# the downloaded file is float32 data format
# for plotting, we need uint8 data format
# plotting the image with full extent
pygge.easy_plot(thisfile, ax=ax1, bands=[4,2,1], percentiles=[0,99])
WARNING:rasterio._env:CPLE_AppDefined in s2.tif: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples. WARNING:rasterio._env:CPLE_AppDefined in TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.
The coordinate reference system (CRS) of the downloaded image composite is not in the UK national map projection. It is reprojected.
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# printing the EPSG code of our shapefile into which we want to reproject the TCI images
print("Reprojecting image to EPSG projection ", epsg)
# making a file name for our new file
warpfile = thisfile.split(sep='.')[0] + '_warped.tif'
print("We are in this directory: ")
!pwd
print("Input file: ", thisfile)
print("Output file: ", warpfile)
# calling the easy_warp function
tmp = pygge.easy_warp(thisfile, warpfile, epsg)
Reprojecting image to EPSG projection 4326 We are in this directory: /content/work/download
WARNING:rasterio._env:CPLE_AppDefined in s2.tif: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples. WARNING:rasterio._env:CPLE_AppDefined in TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.
Input file: s2.tif Output file: s2_warped.tif Creating warped file:s2_warped.tif
to see the locations of our polygons on top of our image composite,We do that with the Geopandas library.
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# creating a figure with subplots
fig, ax = plt.subplots(2,1, figsize=(10,16))
fig.patch.set_facecolor('white')
# plotting the image with full extent
pygge.easy_plot(warpfile, ax=ax[0], percentiles=[0,98], bands=[3,2,1],
shapefile=shapefile, fillcolor="none", linecolor="yellow",
title="LARACHE AFTER FIRE")
# plotting the image with full extent
pygge.easy_plot(warpfile, ax=ax[1], percentiles=[0,98], bands=[4,2,1],
shapefile=shapefile, fillcolor="none", linecolor="yellow",
title="LARACHE AFTER FIRE")
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# reading in the warped before fire downloaded image of larache and checking its cordinate reference system
larache_before = rasterio.open(join(rootdir + "/larache_before/before_fire_warped.tif"))
src1=larache_before
d1=src1.crs.to_wkt()
print(type(d1))
pprint(d1)
<class 'str'>
('GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS '
'84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]')
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# reading in the warped after fire downloaded image of larache and checking its cordinate system
larache_after = rasterio.open(join(rootdir + "/larache_after/after_fire_warped.tif"))
src2=larache_after
d2=src2.crs.to_wkt()
print(type(d2))
pprint(d2)
<class 'str'>
('GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS '
'84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]')
# checking if the before and after image data are the same shape?
larache_before.shape == larache_after.shape
True
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
swir_before = larache_before.read(5) # band in position 5 in our stacked image i.e SWIR
nir_before = larache_before.read(4) # band in position 4 in our stacked image i.e NIR
# calculation of NBR before the fire as flaoting point array
nbr_before = (nir_before.astype(float) - swir_before.astype(float)) / (nir_before.astype(float) + swir_before.astype(float))
# Ignoring division by zero
np.seterr(divide='ignore', invalid='ignore')
# printing some image statistics, ignoring missing values (nan)
print("minimum NBR_before = ", np.nanmin(nbr_before))
print("mean NBR_before = ", np.nanmean(nbr_before))
print("maximum NBR_before = ", np.nanmax(nbr_before))
print("standard deviation = ", np.nanstd(nbr_before))
minimum NBR_before = -0.14298390392538413 mean NBR_before = 0.10239791290152904 maximum NBR_before = 0.44895833342572644 standard deviation = 0.0735262823032699
plt.imshow(nbr_before, cmap='PiYG') # colours to use include Purple and Green
plt.colorbar()
plt.title('NBR of Larache before fire') # title of plot
plt.show() # plot
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
swir_after = larache_after.read(5) # band in position 5 in our stacked image i.e SWIR
nir_after = larache_after.read(4) # band in position 4 in our stacked image i.e NIR
# calculation of NBR after the fire as flaoting point array
nbr_after = (nir_after.astype(float) - swir_after.astype(float)) / (nir_after.astype(float) + swir_after.astype(float))
# Ignoring division by zero
np.seterr(divide='ignore', invalid='ignore')
# printing some image statistics, ignoring missing values (nan)
print("minimum NBR_after = ", np.nanmin(nbr_after))
print("mean NBR_after = ", np.nanmean(nbr_after))
print("maximum NBR_after = ", np.nanmax(nbr_after))
print("standard deviation = ", np.nanstd(nbr_after))
minimum NBR_after = -0.3224225354386313 mean NBR_after = 0.049028595859026615 maximum NBR_after = 0.5460599474035327 standard deviation = 0.11624440093712922
plt.imshow(nbr_after, cmap='PiYG') # Purple and Green selected as colours in plot
plt.colorbar()
plt.title('NBR of Larache after fire') # title of plot
plt.show() # display plot or map
# finding the difference between the NBR before and after the fire
dNBR = nbr_before - nbr_after
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# printing some image statistics, ignoring missing values (nan)
print("minimum dNBR = ", np.nanmin(dNBR))
print("mean dNBR = ", np.nanmean(dNBR))
print("maximum dNBR = ", np.nanmax(dNBR))
print("standard deviation = ", np.nanstd(dNBR))
minimum dNBR = -0.363571844669598 mean dNBR = 0.05336931704250243 maximum dNBR = 0.5638484811508442 standard deviation = 0.10752007563079775
# Defining the dNBR values and the corresponding severity levels
dNBR
levels = [-0.5, -0.25, -0.1, 0.1, 0.27,0.44,0.66]
labels = ['Enhanced regrowth,high','Enhanced regrowth,low','Unburned', 'Low Severity', 'Moderate-low Severity', 'Moderate-high Severity', 'High Severity']
# Creating a custom color map with 5 colors
cmap = colors.ListedColormap(['purple','cyan','green', 'yellow', 'blue', 'darkred', 'red'])
# Setting the boundaries and norm of the color map
bounds = [-0.5, -0.25, -0.1, 0.1, 0.27,0.44,0.66]
norm = colors.BoundaryNorm(bounds, cmap.N)
# Creating a figure and axis object
fig, ax = plt.subplots(figsize=(6, 6))
# Plotting the dNBR using the custom color map and norm
im = ax.imshow(dNBR, cmap=cmap, norm=norm)
# Adding a color bar with custom tick labels
cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, boundaries=bounds, ticks=levels)
cbar.ax.set_yticklabels(labels)
# Setting the title and axis labels
ax.set_title('dNBR Map with Severity Levels')
plt.show()
# Define the area of each pixel in hectares assuming a spatial resolution of 20 meters
pixel_area = 0.0004
# Calculate the number of pixels in each class or label
regrowth_high_pixels = np.sum((dNBR >= -0.5) & (dNBR < -0.251))
regrowth_low_pixels = np.sum((dNBR >= -0.25) & (dNBR < -0.101))
unburned_pixels = np.sum((dNBR >= -0.100) & (dNBR < 0.99))
low_severity_pixels = np.sum((dNBR >= 0.100) & (dNBR < 0.269))
mod_low_severity_pixels = np.sum((dNBR >= 0.270) & (dNBR < 0.439))
mod_high_severity_pixels = np.sum((dNBR >= 0.44) & (dNBR < 0.659))
high_severity_pixels = np.sum((dNBR >= 0.660) & (dNBR < 1.300))
# Calculate the area of each class or label in hectares
regrowth_high_area = regrowth_high_pixels * pixel_area
regrowth_low_area = regrowth_low_pixels * pixel_area
unburned_area = unburned_pixels * pixel_area
low_severity_area = low_severity_pixels * pixel_area
mod_low_severity_area = mod_low_severity_pixels * pixel_area
mod_high_severity_area = mod_high_severity_pixels * pixel_area
high_severity_area = high_severity_pixels * pixel_area
# Print the area of each class or label
print('Enhanced regrowth,high area:', regrowth_high_area, 'hectares')
print('Enhanced regrowth, low area:', regrowth_low_area, 'hectares')
print('Unburned area:', unburned_area, 'hectares')
print('Low severity area:', low_severity_area, 'hectares')
print('Moderate-low severity area:', mod_low_severity_area, 'hectares')
print('Moderate-high severity area:', mod_high_severity_area, 'hectares')
print('High severity area:', high_severity_area, 'hectares')
Enhanced regrowth,high area: 0.0036000000000000003 hectares Enhanced regrowth, low area: 0.0664 hectares Unburned area: 74.6992 hectares Low severity area: 2.0148 hectares Moderate-low severity area: 4.2816 hectares Moderate-high severity area: 1.4420000000000002 hectares High severity area: 0.0 hectares
# Opening the TIF file and the shapefile of burned polygon
with rasterio.open('/content/drive/MyDrive/satelliteCW1/larache_after/after_fire_warped.tif') as src:
shapes = gpd.read_file('/content/drive/MyDrive/satelliteCW1/burned_layers/POLYGON.shp')
# Masking the TIF file using the shapefile
masked_data, _ = rasterio.mask.mask(src, shapes.geometry, crop=True)
# Calculating the NBR2_burned
# the NIR band has index 3 and SWIR band has index 4
NBR2_burned = (masked_data[3] - masked_data[4]) / (masked_data[3] + masked_data[4])
# Ignoring division by zero
np.seterr(divide='ignore', invalid='ignore')
# print average NBR2 for burned poloygon, ignoring missing values (nan)
print("average NBR2_burned = ", np.nanmean(NBR2_burned))
average NBR2_burned = -0.13870606
# Opening the TIF file and the shapefile of burned polygon
with rasterio.open('/content/drive/MyDrive/satelliteCW1/larache_before/before_fire_warped.tif') as src:
shapes = gpd.read_file('/content/drive/MyDrive/satelliteCW1/burned_layers/POLYGON.shp')
# Masking the TIF file using the shapefile
masked_data, _ = rasterio.mask.mask(src, shapes.geometry, crop=True)
# Calculating the NBR1_burned
# the NIR band has index 3 and SWIR band has index 4
NBR1_burned = (masked_data[3] - masked_data[4]) / (masked_data[3] + masked_data[4])
# Ignoring division by zero
np.seterr(divide='ignore', invalid='ignore')
# printing average NBR1 for burned polygon, ignoring missing values (nan)
print("average NBR1_burned = ", np.nanmean(NBR1_burned))
average NBR1_burned = 0.1715795
CALCULATION OF dNBR FOR BURNED POLYGON
# finding the difference between the NBR before and after the fire
dNBR_burned = NBR1_burned - NBR2_burned
# printing average dNBR for burned polygon, ignoring missing values (nan)
print("average dNBR_burned = ", np.nanmean(dNBR_burned))
average dNBR_burned = 0.3102856
# Opening the TIF file and the shapefile of unburned polygon
with rasterio.open('/content/drive/MyDrive/satelliteCW1/larache_before/before_fire_warped.tif') as src:
shapes = gpd.read_file('/content/drive/MyDrive/satelliteCW1/unburned_layers/POLYGON.shp')
# Masking the TIF file using the shapefile
masked_data, _ = rasterio.mask.mask(src, shapes.geometry, crop=True)
# Calculating the NBR1_unburned
# the NIR band has index 3 and SWIR band has index 4
NBR1_unburned = (masked_data[3] - masked_data[4]) / (masked_data[3] + masked_data[4])
# Ignoring division by zero
np.seterr(divide='ignore', invalid='ignore')
# printing average NBR1 for unburned poloygon, ignoring missing values (nan)
print("average NBR1_unburned = ", np.nanmean(NBR1_unburned))
average NBR1_unburned = 0.17309634
# Opening the TIF file and the shapefile of unburned polygon
with rasterio.open('/content/drive/MyDrive/satelliteCW1/larache_after/after_fire_warped.tif') as src:
shapes = gpd.read_file('/content/drive/MyDrive/satelliteCW1/unburned_layers/POLYGON.shp')
# Masking the TIF file using the shapefile
masked_data, _ = rasterio.mask.mask(src, shapes.geometry, crop=True)
# Calculating the NBR2_unburned
# the NIR band has index 3 and SWIR band has index 4
NBR2_unburned = (masked_data[3] - masked_data[4]) / (masked_data[3] + masked_data[4])
# Ignoring division by zero
np.seterr(divide='ignore', invalid='ignore')
# printing average NBR2 for unburned poloygon, ignoring missing values (nan)
print("average NBR2_unburned = ", np.nanmean(NBR2_unburned))
average NBR2_unburned = 0.17163754
# finding the difference between the NBR before and after the fire
dNBR_unburned = NBR1_unburned - NBR2_unburned
# printing average dNBR for burned polygon, ignoring missing values (nan)
print("average dNBR_unburned = ", np.nanmean(dNBR_unburned))
average dNBR_unburned = 0.0014588113
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
r_before = larache_before.read(3) # band 3 in our stacked image is Red
nir_before = larache_before.read(4) # band 4 in our stacked image is NIR
# calculation of NDVI before the fire as flaoting point array
NDVI_before = (nir_before.astype(float) - r_before.astype(float)) / (nir_before.astype(float) + r_before.astype(float))
# Ignore division by zero
np.seterr(divide='ignore', invalid='ignore')
# print some image statistics, ignoring missing values (nan)
print("minimum NDVI_before = ", np.nanmin(NDVI_before))
print("mean NDVI_before = ", np.nanmean(NDVI_before))
print("maximum NDVI_before = ", np.nanmax(NDVI_before))
print("standard deviation = ", np.nanstd(NDVI_before))
minimum NDVI_before = -0.0031236251459653844 mean NDVI_before = 0.16960801456931313 maximum NDVI_before = 0.44798869776436046 standard deviation = 0.053921411038247144
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
r_after = larache_after.read(3) # band 3 in our stacked image is Red
nir_after = larache_after.read(4) # band 4 in our stacked image is NIR
# calculation of NDVI after the fire as flaoting point array
NDVI_after = (nir_after.astype(float) - r_after.astype(float)) / (nir_after.astype(float) + r_after.astype(float))
# Ignoring division by zero
np.seterr(divide='ignore', invalid='ignore')
# print some image statistics, ignoring missing values (nan)
print("minimum NDVI_after = ", np.nanmin(NDVI_after))
print("mean NDVI_after = ", np.nanmean(NDVI_after))
print("maximum NDVI_after = ", np.nanmax(NDVI_after))
print("standard deviation = ", np.nanstd(NDVI_after))
minimum NDVI_after = -0.06513785137467486 mean NDVI_after = 0.17277971166410844 maximum NDVI_after = 0.5750123655099424 standard deviation = 0.07998356629775627
plt.imshow(NDVI_before, cmap='RdYlGn') # selection of Red to Green colours
plt.colorbar() # colour bar
plt.title('NDVI of Larache before fire') # map title
plt.show() # showing plot
# Defining the pre-fire NDVI levels
NDVI_before
levels = [-1, 0, 0.2, 0.5]
labels = ['barren or non-vegetated','sparse or stressed vegetation','moderate or healthy vegetation', 'dense or very healthy vegetation']
# Creating a custom color map with 5 colors
cmap = colors.ListedColormap(['orange', 'yellow', 'lightgreen', 'darkgreen'])
# Setting the boundaries and norm of the color map
bounds = [-1, 0, 0.2, 0.5]
norm = colors.BoundaryNorm(bounds, cmap.N)
# Creating a figure and axis object
fig, ax = plt.subplots(figsize=(6, 6))
# Plotting the pre-fire NDVI using the custom color map and norm
im = ax.imshow(NDVI_before, cmap=cmap, norm=norm)
# Adding a color bar with custom tick labels
cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, boundaries=bounds, ticks=levels)
cbar.ax.set_yticklabels(labels)
# Setting the title and axis labels
ax.set_title('Pre-fire NDVI Map of Larache ')
plt.show()
MAPS OF NDVI AFTER THE FIRE
plt.imshow(NDVI_after, cmap='RdYlGn') # Red to Green colour selected
plt.colorbar()
plt.title('NDVI of Larache after fire') # title
plt.show()
# Defining the post-fire NDVI levels
NDVI_after
levels = [-1, 0, 0.2, 0.5]
labels = ['barren or non-vegetated','sparse or stressed vegetation','moderate or healthy vegetation', 'dense or very healthy vegetation']
# Creating a custom color map with 5 colors
cmap = colors.ListedColormap(['orange', 'yellow', 'lightgreen', 'darkgreen'])
# Setting the boundaries and norm of the color map
bounds = [-1, 0, 0.2, 0.5]
norm = colors.BoundaryNorm(bounds, cmap.N)
# Creating a figure and axis object
fig, ax = plt.subplots(figsize=(6, 6))
# Plotting the post-fire NDVI using the custom color map and norm
im = ax.imshow(NDVI_after, cmap=cmap, norm=norm)
# Adding a color bar with custom tick labels
cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, boundaries=bounds, ticks=levels)
cbar.ax.set_yticklabels(labels)
# Setting the title and axis labels
ax.set_title('Post-fire NDVI Map of Larache ')
plt.show()
# finding the difference between the NDVI before and after the fire
dNDVI = NDVI_before - NDVI_after
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
# printing some image statistics, ignoring missing values (nan)
print("minimum dNDVI = ", np.nanmin(dNDVI))
print("mean dNDVI = ", np.nanmean(dNDVI))
print("maximum dNDVI = ", np.nanmax(dNDVI))
print("standard deviation = ", np.nanstd(dNDVI))
minimum dNDVI = -0.3098725340997325 mean dNDVI = -0.0031716970947953077 maximum dNDVI = 0.33930637628747506 standard deviation = 0.06478898685069448
# Defining the dNDVI values and the corresponding severity levels
dNDVI
levels = [-0.5, -0.25, -0.1, 0.1, 0.27,0.44,0.66]
labels = ['Enhanced regrowth,high','Enhanced regrowth,low','Unburned', 'Low Severity', 'Moderate-low Severity', 'Moderate-high Severity', 'High Severity']
# Creating a custom color map with 5 colors
cmap = colors.ListedColormap(['purple','cyan','green', 'yellow', 'blue', 'darkred', 'red'])
# Setting the boundaries and norm of the color map
bounds = [-0.5, -0.25, -0.1, 0.1, 0.27,0.44,0.66]
norm = colors.BoundaryNorm(bounds, cmap.N)
# Creating a figure and axis object
fig, ax = plt.subplots(figsize=(6, 6))
# Plotting the dNDVI using the custom color map and norm
im = ax.imshow(dNDVI, cmap=cmap, norm=norm)
# Adding a color bar with custom tick labels
cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, boundaries=bounds, ticks=levels)
cbar.ax.set_yticklabels(labels)
# Setting the title and axis labels
ax.set_title('dNDVI Map with Severity Levels')
plt.show()
# Define the area of each pixel in hectares assuming a spatial resolution of 20 meters
pixel_area = 0.0004
# Calculate the number of pixels in each class or label
regrowth_high_pixels = np.sum((dNDVI >= -0.5) & (dNDVI < -0.251))
regrowth_low_pixels = np.sum((dNDVI >= -0.25) & (dNDVI < -0.101))
unburned_pixels = np.sum((dNDVI >= -0.100) & (dNDVI < 0.99))
low_severity_pixels = np.sum((dNDVI >= 0.100) & (dNDVI < 0.269))
mod_low_severity_pixels = np.sum((dNDVI >= 0.270) & (dNDVI < 0.439))
mod_high_severity_pixels = np.sum((dNDVI >= 0.44) & (dNDVI < 0.659))
high_severity_pixels = np.sum((dNDVI >= 0.660) & (dNDVI < 1.300))
# Calculate the area of each class or label in hectares
regrowth_high_area = regrowth_high_pixels * pixel_area
regrowth_low_area = regrowth_low_pixels * pixel_area
unburned_area = unburned_pixels * pixel_area
low_severity_area = low_severity_pixels * pixel_area
mod_low_severity_area = mod_low_severity_pixels * pixel_area
mod_high_severity_area = mod_high_severity_pixels * pixel_area
high_severity_area = high_severity_pixels * pixel_area
# Print the area of each class or label
print('Enhanced regrowth,high area:', regrowth_high_area, 'hectares')
print('Enhanced regrowth, low area:', regrowth_low_area, 'hectares')
print('Unburned area:', unburned_area, 'hectares')
print('Low severity area:', low_severity_area, 'hectares')
print('Moderate-low severity area:', mod_low_severity_area, 'hectares')
print('Moderate-high severity area:', mod_high_severity_area, 'hectares')
print('High severity area:', high_severity_area, 'hectares')
Enhanced regrowth,high area: 0.004 hectares Enhanced regrowth, low area: 0.2932 hectares Unburned area: 74.45960000000001 hectares Low severity area: 6.163600000000001 hectares Moderate-low severity area: 0.15760000000000002 hectares Moderate-high severity area: 0.0 hectares High severity area: 0.0 hectares
import numpy as np
from scipy.stats import linregress
import matplotlib.pyplot as plt
# Assuming dNBR and dNDVI are 2D arrays with the same shape
dNBR_flat = -dNBR.flatten()
dNDVI_flat = dNDVI.flatten()
# Calculate the linear regression line
slope, intercept, r_value, p_value, std_err = linregress(dNBR_flat, dNDVI_flat)
line = slope * dNBR_flat + intercept
# Create a scatter plot
fig, ax = plt.subplots()
ax.scatter(dNBR_flat, dNDVI_flat, alpha=0.5)
# Add the linear regression line to the scatter plot with the equation and correlation coefficient
eqn_label = f'dNDVI = {slope:.2f}(dNBR) + {intercept:.2f}'
corr_label = f'r = {r_value:.2f}'
ax.plot(dNBR_flat, line, color='red', label=f'{eqn_label}\n{corr_label}')
# Set the x-axis and y-axis labels
ax.set_xlabel('dNBR')
ax.set_ylabel('dNDVI')
# Set the title of the plot
ax.set_title('Plot of dNDVI versus dNBR')
# Add the legend to the plot
ax.legend()
plt.show()
'''
--------------------------------------------------------------
The following block of code originates and is modified from:
Balzter, H. (2023)
Materials for GY7709 masters computer classes.
https://uniofleicester-my.sharepoint.com/personal/hb91_leicester_ac_uk/_layouts/15
/onedrive.aspx?id=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%2FGY7709%5FSatellite%5FData%5FAnalysis%
5Fin%5FPython%2F2022%2D23%2FGY7709%2D2022%2D23%2Ezip&parent=%2Fpersonal%2Fhb91%5Fleicester%5Fac%5Fuk%2FDocuments%2FDesktop%
2FGY7709%5FSatellite%5FData%5FAnalysis%5Fin%5FPython%2F2022%2D23&ga=1
Downloaded 1 February 2023
--------------------------------------------------------------
'''
#To save a notebook as an html file, i add a new code cell at the end of the notebook with the following contents:
!pip install -U notebook-as-pdf
!pyppeteer-install
!jupyter nbconvert /content/drive/MyDrive/satelliteCW1/229010645_GY7709_CW2.ipynb --to html
Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Collecting notebook-as-pdf
Downloading notebook_as_pdf-0.5.0-py3-none-any.whl (6.5 kB)
Requirement already satisfied: nbconvert in /usr/local/lib/python3.10/dist-packages (from notebook-as-pdf) (6.5.4)
Collecting pyppeteer (from notebook-as-pdf)
Downloading pyppeteer-1.0.2-py3-none-any.whl (83 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 83.4/83.4 kB 3.7 MB/s eta 0:00:00
Collecting PyPDF2 (from notebook-as-pdf)
Downloading pypdf2-3.0.1-py3-none-any.whl (232 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 232.6/232.6 kB 11.9 MB/s eta 0:00:00
Requirement already satisfied: lxml in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (4.9.2)
Requirement already satisfied: beautifulsoup4 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (4.11.2)
Requirement already satisfied: bleach in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (6.0.0)
Requirement already satisfied: defusedxml in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (0.7.1)
Requirement already satisfied: entrypoints>=0.2.2 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (0.4)
Requirement already satisfied: jinja2>=3.0 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (3.1.2)
Requirement already satisfied: jupyter-core>=4.7 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (5.3.0)
Requirement already satisfied: jupyterlab-pygments in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (0.2.2)
Requirement already satisfied: MarkupSafe>=2.0 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (2.1.2)
Requirement already satisfied: mistune<2,>=0.8.1 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (0.8.4)
Requirement already satisfied: nbclient>=0.5.0 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (0.7.4)
Requirement already satisfied: nbformat>=5.1 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (5.8.0)
Requirement already satisfied: packaging in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (23.1)
Requirement already satisfied: pandocfilters>=1.4.1 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (1.5.0)
Requirement already satisfied: pygments>=2.4.1 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (2.14.0)
Requirement already satisfied: tinycss2 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (1.2.1)
Requirement already satisfied: traitlets>=5.0 in /usr/local/lib/python3.10/dist-packages (from nbconvert->notebook-as-pdf) (5.7.1)
Requirement already satisfied: appdirs<2.0.0,>=1.4.3 in /usr/local/lib/python3.10/dist-packages (from pyppeteer->notebook-as-pdf) (1.4.4)
Requirement already satisfied: certifi>=2021 in /usr/local/lib/python3.10/dist-packages (from pyppeteer->notebook-as-pdf) (2022.12.7)
Collecting importlib-metadata>=1.4 (from pyppeteer->notebook-as-pdf)
Downloading importlib_metadata-6.6.0-py3-none-any.whl (22 kB)
Collecting pyee<9.0.0,>=8.1.0 (from pyppeteer->notebook-as-pdf)
Downloading pyee-8.2.2-py2.py3-none-any.whl (12 kB)
Requirement already satisfied: tqdm<5.0.0,>=4.42.1 in /usr/local/lib/python3.10/dist-packages (from pyppeteer->notebook-as-pdf) (4.65.0)
Requirement already satisfied: urllib3<2.0.0,>=1.25.8 in /usr/local/lib/python3.10/dist-packages (from pyppeteer->notebook-as-pdf) (1.26.15)
Collecting websockets<11.0,>=10.0 (from pyppeteer->notebook-as-pdf)
Downloading websockets-10.4-cp310-cp310-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl (106 kB)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 106.8/106.8 kB 14.4 MB/s eta 0:00:00
Requirement already satisfied: zipp>=0.5 in /usr/local/lib/python3.10/dist-packages (from importlib-metadata>=1.4->pyppeteer->notebook-as-pdf) (3.15.0)
Requirement already satisfied: platformdirs>=2.5 in /usr/local/lib/python3.10/dist-packages (from jupyter-core>=4.7->nbconvert->notebook-as-pdf) (3.3.0)
Requirement already satisfied: jupyter-client>=6.1.12 in /usr/local/lib/python3.10/dist-packages (from nbclient>=0.5.0->nbconvert->notebook-as-pdf) (6.1.12)
Requirement already satisfied: fastjsonschema in /usr/local/lib/python3.10/dist-packages (from nbformat>=5.1->nbconvert->notebook-as-pdf) (2.16.3)
Requirement already satisfied: jsonschema>=2.6 in /usr/local/lib/python3.10/dist-packages (from nbformat>=5.1->nbconvert->notebook-as-pdf) (4.3.3)
Requirement already satisfied: soupsieve>1.2 in /usr/local/lib/python3.10/dist-packages (from beautifulsoup4->nbconvert->notebook-as-pdf) (2.4.1)
Requirement already satisfied: six>=1.9.0 in /usr/local/lib/python3.10/dist-packages (from bleach->nbconvert->notebook-as-pdf) (1.16.0)
Requirement already satisfied: webencodings in /usr/local/lib/python3.10/dist-packages (from bleach->nbconvert->notebook-as-pdf) (0.5.1)
Requirement already satisfied: attrs>=17.4.0 in /usr/local/lib/python3.10/dist-packages (from jsonschema>=2.6->nbformat>=5.1->nbconvert->notebook-as-pdf) (23.1.0)
Requirement already satisfied: pyrsistent!=0.17.0,!=0.17.1,!=0.17.2,>=0.14.0 in /usr/local/lib/python3.10/dist-packages (from jsonschema>=2.6->nbformat>=5.1->nbconvert->notebook-as-pdf) (0.19.3)
Requirement already satisfied: pyzmq>=13 in /usr/local/lib/python3.10/dist-packages (from jupyter-client>=6.1.12->nbclient>=0.5.0->nbconvert->notebook-as-pdf) (23.2.1)
Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.10/dist-packages (from jupyter-client>=6.1.12->nbclient>=0.5.0->nbconvert->notebook-as-pdf) (2.8.2)
Requirement already satisfied: tornado>=4.1 in /usr/local/lib/python3.10/dist-packages (from jupyter-client>=6.1.12->nbclient>=0.5.0->nbconvert->notebook-as-pdf) (6.3.1)
Installing collected packages: pyee, websockets, PyPDF2, importlib-metadata, pyppeteer, notebook-as-pdf
Successfully installed PyPDF2-3.0.1 importlib-metadata-6.6.0 notebook-as-pdf-0.5.0 pyee-8.2.2 pyppeteer-1.0.2 websockets-10.4
[INFO] Starting Chromium download.
100% 109M/109M [00:00<00:00, 234Mb/s]
[INFO] Beginning extraction
[INFO] Chromium extracted to: /root/.local/share/pyppeteer/local-chromium/588429
[NbConvertApp] Converting notebook /content/drive/MyDrive/satelliteCW1/229010645_GY7709_CW2.ipynb to html
[NbConvertApp] Writing 6741948 bytes to /content/drive/MyDrive/satelliteCW1/229010645_GY7709_CW2.html